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The Born effective charge, Z* , that describes the polarization created by collective atomic dis- 
placements, can be computed from first-principles by different techniques. Its band-by-band decom- 
position appeared recently as a powerful tool for the microscopic characterisation of the bonding in 
solids. We describe the connections between the different expressions used for the computation of 
Z* , and analyze the possible associated band-by-band decompositions. We show that unlike for the 
full Z* , the different band-by-band values are not equal, and emphasize that one of them has an 



OO 

I interesting physical meaning. 



I. INTRODUCTION 



The Born effective charge (Z*) describes the macroscopic polarization induced, in a crystalline insulating solid, 
by the collective displacements of nuclei belonging to a given sublattice. In the lattice dynamics study of insulating 
crystals, this tensor is usually considered as a fundamental quantity jjj, because it governs the amplitude of the 
long-range Coulomb interaction between nuclei, and the splitting between longitudinal (LO) and transverse (TO) 
^vq optic phonon modes. 

t-H | In simple materials, like A N B S ~ N binary crystals 0, in which phonon eigenvectors are imposed by symmetry, 
infra-red measurement of the splitting between LO and TO modes allows an accurate estimation of \Z*\ 2 /eoo and 
offers an unambiguous way to extract the amplitude of Z* (its sign remains undefined) from the experiment. However, 
in more complex materials like ABO3 compounds, LO and TO mode eigenvectors are not necessarily equivalent and 
the determination of Z* from the experimental data is not so straightforward : the value can only be approximated 
within some "realistic" hypothesis j3|. For such compounds, the development of theoretical methods giving directly 
' access to Z* was therefore particularly interesting. 

Being defined in term of response to atomic displacements, the Born effective charge is a dynamical concept that 
cannot be assimilated to conventional static charges 0. In particular, the amplitude of Z* cannot be estimated from 
the inspection of the ground-state electronic density alone : It is also dependent on the electronic current flowing 
through the crystal when nuclei are displaced. In some materials, this current may be anomalously large and the 
O ■ amplitude of Z* can deviate substantially from that of the nominal ionic charges |5|-|l(J . 

During the early seventies, the physics of the Born effective charges was already widely discussed within various 
semi-empirical models j|,[ll]-|l4] . Since that time, different theoretical advances have been performed toward their 
hrst-principles determination. A hrst powerful and systematic procedure was introduced by Baroni, Giannozzi and 
Testa |jl5| , who suggested to determine Z* from a linear response formalism grounded on a Sternheimer equation. A 
different algorithm, based on a variational principle, was then reported by Gonze, Allan and Teter |]l6[| , yielding a new 
expression for Z* . Thanks to recent progress in the theory of the macroscopic polarization, Z* is now also directly 
accessible from finite difference of polarization [jlTj . If the first two methods were exclusively implemented within the 
density functional formalism (DFT), the last one also allowed calculations of changes in polarization within different 
other one-electron schemes (Hartree-Fock method model GW approximations to many-body theory p^BO] , 

Harrison tight-binding model El]]) and the Hubbard tight-binding model p2| . 

The Born effective charges are now routinely computed within DFT and accurate prediction have been reported for a 
large variety of materials. Without being exhaustive, these calculations concern monoatomic crystals as selenium [ p3| , 
various I- VII |plp|,I I-VI prl|27| III-V [p8|j29| IV-IV [0Jl_or IV- VI || semiconductors as well as a large diversity 



of oxydes (AO |3j,|3j, A0 2 |6[|Q, A0 3 |0), AB0 3 @ : |7|,|37|1) and even more exotic materials as Al 2 Ru @. The 
previous results exclusively concern crystalline materials but recently a calculation of Z* in a model of amorphous 
silica has also been reported J3S| ]. 

Going beyond the bare determination of Z*, the first-principles approach was also offering a new opportunity to 
investigate the microscopic mechanisms monitoring its amplitude. In some recent studies [ p9|j4C| , pj| , |35| , the decompo- 
sition of individual contributions from separate groups of occupied bands to Z* appeared particularly useful : The 
amplitude of Z* can be related to dynamical changes of hybridizations and the band-by-band decomposition allowed 
to identify the orbitals involved in such a mechanism. 
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In the present paper, we aim at presenting the links between the theoretical frameworks used nowadays for the 
computation of Born effective charges, and deduce from this comparison the correct way to develop a band-by- 
band analysis. We note that unitary transformations among occupied states, called "gauge tranformations" are 
present at different levels in the formalism. These unitary transforms need to be treated carefully in order to obtain 
physical results. Also, the computation of a band-by-band decomposition might be conceived from responses to an 
homogeneous electric field or to a collective ionic displacement. However, some care is needed here in order to yield 
physically meaningful quantities. The final expression that we propose is directly linked to displacements of the center 
of gravity of Wannier functions. 

In a first section, we make a brief overview of the different procedures that, historically, have been considered to 
determine Z* from first-principles. In section II, we detail how Z* can be computed within DFT in the Berry phase 
approach and the linear response formalism (either from the response to an homogeneous electric field or a collective 
ionic displacement). We take into account the possible gauge transformations. Starting from the general expressions 
of Z* , we then describe in Section III how contributions of isolated sets of bands can be separated from each others. 
We discuss the physical significance of such band-by-band contributions in terms of Wannier functions. We illustrate 
our discussion with a band- by-band decomposition of the titanium charge in BaTiC>3, using different expressions. 



II. OVERVIEW OF THE DIFFERENT THEORETICAL FORMULATIONS OF Z 



A. Equivalent definitions of Z* 

For periodic systems, the Born effective charge tensor Z* a o of nuclei belonging to the sublattice k is conventionally 
defined as the coefficient of proportionality relating, under the condition of zero macroscopic electric field, the change in 
macroscopic polarization Vp along the direction (3 and the collective nuclear displacements of atoms k along direction 
a, times the unit cell volume £l : 



z* -n dVp 

U 1 Km 



(l) 

£=0 



We note the important condition of vanishing macroscopic electric field, necessary to fix unambiguously the amplitude 
of the change of polarization. Similarly, other dynamical charges may be defined from the change of polarization 
induced by a sublattice displacement when imposing other conditions on the field £. These charges (Callen charge [ ffl| , 
Szigeti charge ]42|| ) can be related to Z* thanks to expressions involving the dielectric tensor (see, for instance, 
Ref. [[§). 

The standard definition of Z* , Eq. (0), emphasizes the response with respect to the collective nuclear displacement. 
However, a thermodynamical equality relates the macroscopic polarization to the derivative of the electric enthalpy 
E with respect to a homogeneous electric field. Similarly another relationship connects the forces on the nuclei to 
the derivative of the electric enthalpy with respect to atomic displacements. Combining these expressions, Z* can be 
formulated, either as a mixed second-order derivative of the electric enthalpy, 

Z* - d ^ (2) 

or as the derivative of the force felt by a nucleus k with respect to an homogeneous effective electric field £p, at zero 
atomic displacements : 



6F-, 



d£ s 



(3) 

T Ka =0 



The three previous definitions - Eqs. (^), (|J), and (^) - are formally equivalent. However, as it is now described, 
each of them can lead to different algorithms for the computation of Z* . 



B. The computation of Z* based on finite differences of the polarization 



From Eq. (|]J), it appears that an easy access to the macroscopic polarization would allow the computation of Z* 
by finite difference. Starting from the equilibrium positions of the atoms and considering a small but finite collective 
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displacement Ar KjCt in direction a of the atoms belonging to the sublattice k, we obtain the following finite-difference 
estimate of Z* from the change of polarization in zero-field : 

Zl a(3 = fi lim £?2- (4) 

In the earliest first-principles computation of Z* , in 1972, Bennett and Maradudin [Q had attempted to deduce Z* 
using this technique, but as pointed out by Martin p5[ , on the basis of an incorrect expression for the polarization, 
yielding boundary-sensitive results. The basic problem was that the change of polarization was assimilated to the 
change of the unit cell dipole, which is ill-defined for delocalized periodic charge distributions. A supercell technique, 
allowing correct estimation of the polarization induced by longitudinal atomic displacements, was then proposed by 
Martin and Kunc Q,^7| to compute the longitudinal charge, giving indirectly access to Z* thanks to the knowledge 
of the optical dielectric constant. 



In 1992, Resta 48 formulated the change of polarization as an integrated macroscopic current. This yielded King- 
Smith and Vandcrbilt [jl7|] to identify in the change of electronic polarization a geometric quantum phase (or Berry 
phase). They proposed a new scheme, useful for a practical calculation of the macroscopic polarization jl9| and from 
which Z* can be deduced by finite difference, from single cell calculations. It is usually referred to as the Berry phase 
approach. 



C. The computation of Z* based on linear responses 

Alternatively, the derivatives can also be computed explicitely using perturbation theory. During the early seventies, 
the linear response of solids to collective nuclear displacements was formulated in term of the inverse dielectric matrix 
e _1 . Accordingly, an expression was proposed for calculating the Born effective charges |50|-|52|. Computations based 
on this formalism were, for example, reported by Resta and Baldereschi [ p6|]53f| . However, one important drawback 
of this procedure was the difficulty to control and to guarantee the charge neutrality, which imposes constraints on 



the off-diagonal elements of e 54 1. Consequently, Vogl |>4j] proposed a method that bypasses the inversion of the 
dielectric function by using directly the self-consistent potential induced by a long- wavelength lattice displacement. 
Unfortunately, at that time, there was no way to determine accurately this potential : it had to be approximated and 
this formulation was only applied to simplified models [5^j5^ ] . 

A solution to this problem was reported in 1987 by Baroni, Giannozzi and Testa |l5| , p8[ ] who proposed, within DFT, 
to compute the total effective potential by solving a self-consistent first-order Sternhcimer equation. It was the first 
"ab initio" powerful and systematic approach yielding accurate value of Z* . Both types of perturbation (the response 
to an homogeneous electric field or to a collective ionic displacement) could be considered in this approach, giving 
access to expressions linked to Eqs. (Q) and (||). A variational formulation of this theory was then reported by Gonze, 



Allan and Teter jl6|,|57 58 1, offering a different algorithm for the calculation of the first-order wavefunctions and a 
new, stationary, expression for Z* a g, linked to Eqs. (|2|). The "Sternheimer" and "variational" formalisms were first 
implemented within DFT when using plane-wave and different kind of pseudopotentials p8P, |5l^ , |6Cfl (usually within 
the LDA, but also within the GGA |6lJ]). Later, LMTO |62j and LAPW versions of the linear response approaches 
have also been proposed. 



D. The computation of Z* based on finite differences of the atomic forces 

In a similar spirit to the Berry phase approach, starting from Eq. (^|), Z* might also a priori be estimated from the 
force linearly induced on the different atoms by a finite homogeneous electric field, when keeping the atomic positions 
frozen : 

The force felt by an atom can be trivially computed using the Hcllmann-Feynman theorem. However, the basic problem 
inherent to this method is that infinite periodic systems (as those considered in practical calculations) cannot sustain 
a finite linear change of potential. Different supercell techniques, that circumvent the problem, have been proposed 
to compute Z* by Kunc and Martin p7[ (from the force in the depolarizing field associated to a longitudinal phonon) 
and by Kunc and Resta |m|,|35| (from the force in an applied sawlike potential). Their calculations demonstrate the 
feasability of the approach. Up to now, it was however only marginally applied. For that reason, in what follows, we 
will concentrate on the linear response formalism and the Berry phase approach. 
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III. ANALYSIS OF THE MATHEMATICAL LINKS BETWEEN THE DIFFERENT FORMULATIONS 



A. Ionic and electronic contributions 



For practical purposes, the Born effective charge can be conveniently decomposed into two contributions : 

%K,aP = Z K 5 a {) + Z% >a p. (6) 

The first term, Z K , is the charge of the nuclei (or pseudo-ion, in case of pseudopotential calculations), and can be 
trivially identified. The second, Z^ a p, is the contribution due to the electrons, in response to the perturbation. In this 
Section we analyse the mathematical links existing between the different expressions that can be used to determine 
this last contribution within the density functional formalism. Similar expressions should be obtained within other 
frameworks as the Hartree-Fock method. 



B. The Berry phase approach 

As previously mentioned, a straightforward approach to the determination of Z* K consists in computing the difference 
of macroscopic polarization between a reference state, and a state where the atoms belonging to the sublattice k have 
been displaced by a small but finite distance Ar K)CB . Combining Eqs. and (||), the electronic contribution to this 
charge can be obtained as : 

AV% 1 

Z°[ a p=n lim ^-S- (7) 

As demonstrated recently by King-Smith and Vanderbilt |l7j| , in periodic systems, the change in electronic polarization 
in zero field can be deduced from the following formula : 

1 OCC. p q 

v t = -^% s ]J u ^w^ ) dk (8) 

where s is the occupation number of states in the valence bands (s = 2 in spin-degenerate system) and u„k is the 
periodic part of the Bloch functions. Taken independently at each wavevector k in the Brillouin zone, the matrix 
elements of the previous equation are ill-defined. Indeed, the phase of the wavefunctions at each wavevector is 
arbitrary, and thus unrelated with the phases at neighbouring wavevectors : The derivative of the wavefunctions with 
respect to their wavevector has a degree of arbitrariness. However, the integral of the right-hand side is a well-defined 
quantity, which has the form of a Berry phase of band n, as discussed by Zak |6^| . 

The King-Smith and Vanderbilt definition is valid only if the wavefunctions fulfill the periodic gauge condition. This 
means that the periodic part of Bloch functions must satisfy 

%k(r) = e lGr u nk+G {r), (9) 

This condition does not fix unambiguously the phase of the wavefunctions at a given k-point (even not at neighbouring 
k-points) but it imposes a constraint between wavefunctions at distant wavevectors. It defines a topology in k-space, 
within which the polarization takes the convenient form of a Berry phase. 

When working within one-electron schemes (DFT, Hartree-Fock, ...), a generalized choice of phase is also present 
at another level. For the ground-state, the Lagrange multiplier method applied to the minimization of the Hohenbcrg 
and Kohn fonctional under orthonormalization conditions on the wavefunctions J67| , gives the following equations : 

occ 

H-^UmSa) = ^ Amn.k |Unk) (10) 

n 

This condition, associated with the minimisation of the Hohenberg and Kohn energy functional, means that the 
wavefunctions u„k must define a Hilbert space (right-hand side of Eq. (|To|)), in which the vector generated by the 
application of the Hamiltonian to one of these wavevectors (left-hand side of Eq. (|l(])) must lie. We observe that a 
unitary transformation between the wavefunctions will leave the Hilbert space invariant, and Eq. (|l0|) will be satisfied 
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provided the matrix of Lagrange multiplier A mjli k is transformed accordingly. In order to build Kohn-Sham band 
structures, the unitary transform is implicitely chosen such as to guarantee 

Amn,k — 0~rnn ^m,k (H) 

in which case 6771, k correspond to the eigenvalues of the Kohn-Sham Hamiltonian and the associated functions Ud,m.k 
are the Kohn-Sham orbitals. This choice is called the diagonal gauge condition. All along this work, it will be 
emphasized by a "d" subscript. 

We note that the periodic gauge condition connect wavefunctions at different k-points, while the diagonal gauge 
condition fixes wavefunctions at a given k-point. The choice defined by Eq. ( pi] ) is not mandatory, and the computation 
of the total energy, the density, or the Berry phase, Eq. (|J), will give the same value independently of the fulfillment 
of Eq. ([n]). If the diagonal gauge is the natural choice for the ground-state wavefunctions, another choice is usually 
preferred for the change in wavefunctions in linear- response calculations (see below). 

From a practical viewpoint, direct evaluation of Eq. (^|) is not trivial in numerical calculations because the wave- 
functions are only computed at a finite number of points in the Brillouin zone, without any phase relationship between 
the eigenvectors. An elegant scheme to deal with this problem was reported in Ref. |l7j]. Also, we note that, associ- 
ated to the fact that a phase is only defined modulo 2n, Eq. (||) only provides the polarization modulo a "quantum" 
(in 3-dimcnsional solids, the quantum is (s.R a /Qo), where R a is a vector of the reciprocal lattice). This quantum 
uncertainty should appear as a limitation of the technique. However, for the purpose of computing Z* from finite- 
differences, the atomic displacements Ar KjQ may always be chosen sufficiently small for the associated change of 
polarization being unambiguously defined. 

The Berry phase technique was successfully applied to various ABO3 compounds HJtJ, giving results equivalent to 
those obtained independently using perturbative techniques |||i8| . A similar agreement was for instance reported for 
alkaline-earth oxides J34|,|35| . 



C. Z* as the first derivative of the polarization 

Instead of approximating Eq. (|l|) from finite differences, it can be chosen to compute it directly. The combination 
of Eqs. (|), (§) and (§) gives : 

1 p (2tt) j ^ J BZ dr K . a dkp dkp dr K . a 

The second expectation value can be worked out : 

, d du nk 

BZ dk/3 OT K . a 

, d ,du nk du nk <9u„ k 

("nklTj } _ (-5i Iq )J« k ( 13 ) 



BZ dkp <9t KiQ dkp dr K ^' 

In the previous expression, the first term of the right-hand side is the gradient of a periodic quantity integrated 
over the Brillouin zone. Within any periodic gauge, Eq. ([)]), its contribution will be zero. Using the time-reversal 
symmetry, we arrive therefore at the final expression : 

Z ^ = ' 2 (^¥ l ? s Lfe 1 w >dk (14) 

The first-derivatives of the wave functions, — u„k and gf^-w n k, required to evaluate this expression, can be 
computed by linear-response techniques either by solving a first-order Sternheimer equation |l5|,^8| or by the direct 
minimization of a variational expression as described in Ref. p^ , |57[ |. 

We note that the choice of gauge will influence the value of the first-derivative of w„k, although the integrated 
quantity Z^ a g must remain independent of this choice (in any periodic gauge). Usually the following choice is 
preferred in linear-response calculations : 



du nil 



Kk) = (15) 

p 
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for m and n labelling occupied states, and A representing either the derivative with respect to the wavevector or to 
atomic displacements. As emphasized by the notation this condition defines what is called the parallel gauge 

and insures that the changes in the occupied wavefunctions are orthogonal to the space of the ground-state occupied 
wavefunctions. This projection on the conduction bands is not reproduced within the diagonal gauge defined by the 
generalization of Eq. ( |il"| ) at the first order of perturbation, as elaborated in Ref. |67j] . 



D. Z* as a mixed second derivatrive of the electric enthalpy 



Considering now Eq. 



Z* appears also as a mixed second derivative of the electric enthalpy. Therefore, following 



the formalism used in Ref. [ p7||58| , Z^ a ^ can be formulated in terms of a stationary expression, involving the first- 
order derivative of the wavefunctions with respect to a collective displacement of atoms of the sublattice k and the 
first-order derivatives of the wavefunctions with respect to an electric field and to their wavevector [|69| : 



J K,a/3 



= 2 



(27T) 



BZ 



E' 



, dUn 



, du r , 



, ,<9umk du 
+(- 1« 



1 f ^dVxcO 

1 



dr K 



d£ 6 



, (r)][^(r)]*dr 

dr Ka dtp 

dn 



+^L,^ F [^(G)]-^(G) 



+\l ^(r,r')[^5-(r)]«^-(r')<irdr' 

1 JClo OT Ka OL(3 

The stationary character of Eq. (fig) is related to the influence of inaccuracies of the derivatives 



(16) 



i)T K 



and 



on the accuracy of Z^ a o : the error in the latest is proportional to the product of the errors in the two wavefunction 
derivatives. By contrast, Eq. ( |l4| ) does not have this property : in that case, the errors on Z* is directly proportional 
to the error on 



jk. As briefly discussed below, Eqs. (|T^) and ( |l4| ) are mathematically equivalent. However, if 
both are estimated with approximate wavefunctions, Eq. ( |l6| ) is more accurate than Eq. (|l4|). 

Eq. ( pT[ ) can be seen to arise from the stationary property of Eq. ([r]). Indeed, in the latter, the error on Z* is 
proportional to the product of the errors on the first-order change in wavefunctions so that if — ti n k was known 

w nk . Putting 



perfectly, a correct estimation of Z^ a0 should be obtained independently of the knowledge of 



dE f 3 



therefore to zero g^w n k and the corresponding density changes in Eq. (|16[), most of the terms cancel out and 
we recover Eq. (|l4|), which evaluated for the exact — « n k must still correspond to a valid expression for Z* . 
Alternatively, the formal equivalence between Eqs. ( |l4| ) and (16) should also be trivially established when introducing 
in Eq. (Rm the first-order Sternheimer equation associated to the atomic displacement perturbation. 



E. Z* as the first derivative of the atomic force 



By the same token as in the preceding Section, we can choose alternatively for — it n k and the associated density 
derivative to vanish in Eq. (flq), and we still obtain a valid expression for Z* : 




This last equation corresponds to the third formulation of Z* in which it appears as the first derivative of the force 
on the atoms k with respect to an electric field (Eq. (||)). Indeed, it is directly connected to the following expression 
of the force, deduced from the Hellmann-Feynman theorem : 
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T^el ll / / I ext,k i . 



(2tt) 3 

+ / [|^(r)][n(r)]*dr (18) 
Jn oT Ka 

Compared to Eqs. and (|l6|), Eq. ([l7|) has the advantage that the computation of the first-order wavefunction 
derivative with respect to the electric field perturbation is the only computationally intensive step needed to deduce 
the full set of effective charges. 

We note however that the implementation of Eq. ( |l7| ) , rather easy within a plane wave - pseudopotential approach, 
is not so straightforward when the basis set is dependent on the atomic positions, as in LAPW methods (additional 
Pulay terms must be introduced). 



IV. A MEANINGFUL BAND-BY-BAND DECOMPOSITION 



A. Problematics 

In the previous Section, we have described different approaches for computing Z^ a g. In order to have a better 
insight on the mechanisms monitoring the amplitude of Z* , it is also useful to separate the partial contribution of 
electrons from the different occupied bands to Z^ a ^. In what follows, we describe how to identify the change of 
polarization induced by the electrons of one specific band and we discuss the meaning of this contribution in terms of 
Wannicr functions. 



B. The displacement of center of gravity of Wannier functions 

Inspired by a previous discussion by Zak p6| , Vanderbilt and King-Smith |70| emphasized that the macroscopic 
electronic polarization acquires a particular meaning when expressed in terms of localized Wannier functions. The 
periodic part of Bloch functions M„k(r) are related to the Wannier functions W n {r) through the following transfor- 
mation : 



XV 



)= e^ k '( r - R ) W n (r - R) (19) 

(20) 



(2tt 

From this definition, we deduce that : 



W n(r) = I e 4k r u„ k (r) dk (21) 

I BZ 



u nk (r) = J= Hfo ~ fy)] e^ k ' (r ~ R) W n (r - R) (22) 



d 

where R runs over all real space lattice vectors. Introducing this result in the equation providing Pp (Eq. ||), we 
obtain : 

V t= £/ rp.\W n (r)\ 2 dr (23) 

° n 

From this equation, the electronic part of the polarization is simply deduced from the position of the center of gravity 
of the electronic charge distribution, as expressed in terms of localized Wannier functions. In other words, for the 
purpose of determining the polarization, "the true quantum mechanical electronic system can be considered as an 
effective classical system of quantized point charges, located at the centers of gravity associated with the occupied 
Wannier functions in each unit cell " [|7Q ], 

We observe that Eqs. (|l^) and (|21|) establish a one-to-one correspondence between u„k and W n . As previously 
emphasized in Section III, when working within the diagonal gauge, Ud.nk becomes identified with the Kohn-Sham 
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orbitals so that the associated W d , n w iU correspond to a single band Wannier function. Within this specific gauge, 
we can therefore isolate P m ,0, the contribution of band m to the [3 component of the polarization, by separating the 
different term in the sum appearing in Eq. (p3|) : 



■pel 



rp.\W d ,m(r)\ 2 dr 



(24) 



If we take now the derivative of the polarization with respect to a collective atomic displacement, Z^ aS can be 
written in terms of Wannier functions as : 



occ 



J n,af3 



rp[( 



dW n (r) 



dr K 



)* W d , n {r) + (W d , n (r)Y 



dW n (r) 



dr K , a 



]dr 



(25) 



As for the polarization, this equation has a simple physical meaning. In response to an atomic displacement, the 
electronic distribution is modified and the electronic contribution to Z* can be identified from the displacement of the 
center of gravity of the occupied Wannier functions. Working within the diagonal gauge at any order of perturbation, 
we will be able to follow the change of single band Wannier functions all along the path of atomic displacements. In 
the previous expression, the contribution of band m to Z e K a g can therefore also be isolated : 



yel ] 
J K,,a(3]m 



dW m (r) 



dr K . 



W d , m (r) + (W d , m (T)y 



dW m (r) 



}dr 



(26) 



This last equation identifies the contribution from band m to the Born effective charge as fio times the change of 
polarization corresponding to the displacement of a point charge s on a distance equal to the displacement of the 
Wannier center of this band. 
Eq. (|26|) can also easily be evaluated in reciprocal space : 



tin 



(2tt) 3 



du r< 



BZ 



dk 







)dk 



(27) 



As Bloch and Wannier functions were related through a band-by-band transformation, the contribution from band m 
to Z* a8 in Eq. ( p7| ) keeps the same clear physical meaning as in Eq. ( |2^ ) : [Z^ l aj3 ] m corresponds to f2 . AV^ 3 — 
£lo ( s.Adp) where A.dp is the displacement in direction fj of the Wannier center of band m induced by the unitary 
displacement of the sublattice of atoms re in direction a. In practical calculations, where each band may be thought as 
a combinaison of well-known orbitals, the displacement of the Wannier center is associated to the admixture of a new 
orbital character to the band and must be attributed to dynamical changes of orbital hybridizations. As illustrated 
in some recent studies |3^^,|l^,^5[ , the decomposition of Z* appears therefore as a powerful tool for the microscopic 
characterisation of the bonding in solids. 

Let us emphasize again that the previous decomposition in terms of a single band is valid only if the diagonal 
gauge was used to define the Kohn-Sham wavefunctions, hence the "d" subscript in Eq. (26) and (j27j). The ground- 
state wavefunctions are conventionally computed within the diagonal gauge. However, in most calculations, the 
first-derivatives of these wavefunctions are computed within the parallel gauge. Within this choice, the change in 
each Bloch function will be a mixing of different Kohn-Sham orbitals when the perturbation is applied so that the 
associated change in functions W n , defined from Eq. (^l|), will correspond to the change of a multi-band Wannier 
function. Evaluating Eq. (|2^) or (|27j ) within such a gauge, we will identify the displacement of a complex of bands 

du n 



but not of a single band. In practice, the first-order derivative of wavefunctions in the diagonal gauge | d can be 
deduced from those in the parallel gauge | p and the ground-state wavefunctions in the diagonal gauge w^nk, by 
adding contributions from the subspace of the occupied bands : 



dUn 



dX 



du 



mk 



d\ 



E 



(Ud : 



I dH |„, 
■tk| gx \ a d,>. 



(28) 



We note that this transformation (Eq. (|2S|)) can present some problems when the denominator vanishes : this happens 
when the valence energies are degenerated. The problem can be partly bypassed by keeping a parallel transport gauge 
within the space of degenerated wavefunctions. Practically, this means that we will only be able to separate the 
contributions of disconnected set of bands. 
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C. Other band- by-band decompositions 



Up to now, we focused on Eq. (fh|). We would like to investigate whether other band-by-band decompositions could 
be obtained from Eqs. (|l6| ) and (p7|). These expressions, unlike Eq. (14), are not written as a simple sum of matrix 
elements, each related with a single band. However, Eq. ([l7|) can still be transformed using a decomposition of the 
density with respect to the bands : 



n(r) 



p OCC 

(o~^3 / <k( r )«nk(r)dk. 

( Z7T ) Jbz _ 



(29) 



It gives 



J n,ai3 



= 2 



^0 



BZ ' 



E, i^cxt.k , dv xc0 du nk 
s {Urik\— — - + l^r^)ak 



dr Ka dr Ka d£ i 







(30) 



for which the following decomposition is obtained, using the diagonal gauge wavefunctions 



a/3\ 



(2tt) 3 



, . dv LtM 

S (Ud,mk|-Q — 

BZ OT Ka 



dr K 



0£ 



(31) 



This expression corresponds to the contribution of the electrons of band m to the force induced on atom k by a 
macroscopic field £p. However, it is not equivalent to Eqs. ( p6| ) or (p7f). Indeed, for a particular band to, the 
difference between matrix elements present in Eq. rt27]) and (ph) is (within a given gauge) : 



d<xt,k 
dr Ka 



dv^rn , du 



mk \ 



dr Ka d£p 
1 

_ 2 



,du 



mk i 



.du 



mk \ 



1 <9tv q 



^c(r,0[^(r)]*^» k 



9r K 



d£p 



(r')rfr dr' 



i / K xc (r,r')[§^(r)]^(r')drdr' 



(32) 



where n TO k(r) is a short notation for w* nk (r)u m k(r). The summation of these differences on all the bands and 
integration on the Brillouin zone gives zero, as expected. However, the band- by-band difference, Eq. (32), does not 
vanish. The quantity defined from Eq. ( |3l|) is therefore independent from that of Eq. ( p6|) and has no specific meaning 
in terms of Wannicr functions. 



D. Numerical comparison of the different decompositions 

The previous theoretical results can now be illustrated on a numerical example. In what follows, we will consider 
the case of barium titanate, a well-known ferroelectric material of the ABO3 perovskite family, presenting non-trivial 
values of Z* . Our calculations have been performed within a planewave-pseudopotential approach. The electronic 
wavefunctions have been expanded in plane- waves up to a kinetic energy cutoff of 35 hartrees. Integrals over the 
Brillouin zone have been replaced by sums on a 6 x 6 x 6 mesh of special k-points. The Born effective charges have 
been obtained by linear response. 

In Table BL we compare different decompositions of the titanium charge (Z^), in the cubic phase. We observe that 
the results obtained within the diagonal and parallel gauge, either from Eq. ( |27| ) or ( |3l| ) are significantly different. 
As demonstrated previously, the identification of meaningful band- by-band contributions require the use of Eq. (|2^), 
when working within the diagonal gauge : the contributions then describe the displacement of the Wannier center of 
each given set of bands, induced in response to the displacement of the Ti atom. A similar decomposition of the Born 
effective charge for the other atoms in the unit cell is reported in Refs. [ 
contributions (i.e. the deviation with respect to the reference ionic values) is discussed in terms of dynamical changes 
of orbital hybridization. 



43|| . In these papers, the origin of anomalous 
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V. CONCLUSIONS 



In this paper, we have recalled the three equivalent definitions of Z* and we have discussed the associated math- 
ematical formulations within the density functional formalism. A particular attention has been placed on the gauge 
freedom associated to unitary transforms within the subspace of the occupied bands. 

Eq. (|l^) expresses Z* as a second derivative of the electric enthalpy. Contrary to the other expressions, it has a 
stationary character and allows the most accurate estimate of Z* when approximate wavefunctions are used. 

Eq. (|l7]) formulates Z* as the force linearly induced on atoms k when a macroscopic electric is applied. It is an 
alternative convenient expression from which the full set of effective charges can be deduced as soon as the derivative 
of the wavefunctions with respect to the electric field perturbation is known. 

Eq. ( |l4| ) and its finite difference expression - Eq.(^) - consider Z* in terms of the macroscopic polarization 
induced by the displacement of the atoms belonging to a given sublattice. Contrary to the other formulations, it 
yields a meaningful band-by-band decomposition, helpful in the characterisation of the bonding in solids. It has been 
demonstrated that the contribution of a particular band m to Z* a/3 , obtained from this expression when working 
within the diagonal gauge, is directly related to the displacement of the Wannier center of this band when a sublattice 
of atoms is displaced. 

Finally, it has been shown explicitely, for BaTiC>3, that the band-by-band decompositions, arising from Eqs. (Q), 
( |l6| ) or (|17|), within the parallel or diagonal gauge, yield significantly different results. 
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TABLE I. Band-by-band decompositions of the Born effective charge of the Ti atom in the cubic phase of BaTi03. In the 
second column are mentioned reference values expected in a purely ionic material; band-by-band contributions presented in 
the three next columns were deduced from first-principles calculations. Only the values obtained from Eq. ( p?] ) and within the 
diagonal gauge can be understood in terms of Wannier functions. 





Reference 
ionic values 


Diagonal gauj 

from Eq. (p7|) 


;e 

from Eq. (|3l|) 


Parallel gauge 
from Eq. (f^) 


core 


+12.00 


+ 12.00 


+12.00 


+12.00 


Ti 3s 


-2.00 


-2.03 


+ 1.56 


-0.36 


Ti 3p 


-6.00 


-6.22 


-9.54 


-5.50 


Ba 5s 


0.00 


+0.05 


-0.36 


0.00 


02s 


0.00 


+0.23 


-1.56 


-0.41 


Ba 5p 


0.00 


+0.36 


+1.47 


+0.10 


O 2p 


0.00 


+2.86 


+3.68 


+1.42 




+4.00 


+7.25 


+7.25 


+7.25 
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